#include "LineEdsCalcTask.h"
#include "mfAddFuncCalcTask.h"
#include "kernel/feObserver.h"
#include <math.h>
#include "Timer.h"
#include "FastMath.h"

//---------------------------------------------------------------------------------------------------------------------------------------------
inline double gCalcPhi(double iX, double iY)
{
	//get phi-coord in polar coords
	if(iX == 0 && iY > 0) return PI/2;
	if(iX == 0 && iY < 0) return 3*PI/2;

	double phi = atan(fabs(iY/iX));
	if(iX >= 0 && iY >= 0) return phi;
	if(iX <= 0 && iY >= 0) return PI - phi;
	if(iX >= 0 && iY <= 0) return 2*PI - phi;	
	if(iX <= 0 && iY <= 0) return phi + PI;

	return phi;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double gCalcR(double iX1, double iY1, double iX2, double iY2)
{
	//calc distance btw 2 points
	iX2 -= iX1;
	iY2 -= iY1;
	return FastMath::Sqrt(iX2*iX2 + iY2*iY2);
}

//---------------------------------------------------------------------------------------------------------------------------------------------
LineEdsCalcTask::LineEdsCalcTask()
{
	NumLayers = 0;
	TimesArray = 0;
	TimesCount = 0;
	DbzDt = new double[1000];
}

//---------------------------------------------------------------------------------------------------------------------------------------------
bool LineEdsCalcTask::Init(ConductiveLayer *iLayers, int iNumLayers, fePolyline<double> *iGenerator, fePolyline<double> *iReceiver)
{
	//create and copy layers
	NumLayers = iNumLayers;
	for(int i=0; i<NumLayers; i++)
		Layers[i] = iLayers[i];
	
	Generator = iGenerator;
	Receiver = iReceiver;

	Dipole.SetMedium(Layers, NumLayers);

	//create triangulation
	RecTris.Init(Receiver);
	RecTris.Triangulate();
	RecTris.Tesselate(1);
	RecTris.FillAreas();
	return true;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
double LineEdsCalcTask::Integrate(double iF[3], double iX[3])
{
/*
	double t1 = iX[2]*iX[2];
	double t2 = t1*iX[2];
	double t14 = iF[0]*t1;
	double t18 = iX[2]*iF[0];
	double t19 = iX[0]*iX[0];
	double t22 = iX[1]*iX[1];
	double t37 = iF[2]*iX[1];
	double t45 = t19*iX[0];
	double t54 = iF[0]*t2-iF[1]*t2+3.0*t1*iF[1]*iX[0]-2.0*t1*iF[2]*iX[1]+2.0*iX[0]*iF[2]*t1
		-4.0*t14*iX[1]+t14*iX[0]-2.0*t18*t19+3.0*t18*t22-iF[2]*t19*iX[2]+2.0*t18*iX[1]*iX[0]+
		3.0*iX[2]*iF[2]*t22-3.0*iF[1]*iX[2]*t19-2.0*iX[2]*iX[0]*t37+4.0*t37*t19-3.0*iF[2]*t22*
		iX[0]+iF[1]*t45-3.0*iF[0]*t22*iX[0]-iF[2]*t45+2.0*iF[0]*iX[1]*t19;
	return t54/(iX[1]-iX[0])/(iX[1]-iX[2])/6.0;
*/


	return 0.5*(iF[0]+iF[1])*(iX[1]-iX[0]);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
bool LineEdsCalcTask::Eds(double iZ, double *iTimes, int iTimesCount, double *oResV)
{
	//check layers
	if(NumLayers < 1 || !Layers || !iTimes || iTimesCount < 1) return false;

	//fill time data
	TimesArray = iTimes;
	TimesCount = iTimesCount;

	ZeroMemory(oResV, sizeof(oResV[0])*iTimesCount);
	
	

	//through lines
	for(int i=0; i<Generator->GetSize(); i++)
	{
		fePolyline<double>::V2 first = Generator->GetPoint(i);
		fePolyline<double>::V2 second = Generator->GetPoint(i+1 >= Generator->GetSize() ? 0 : i+1);
		
		fePolyline<double>::V2 v = second - first;
		double l = FastMath::Sqrt(v.X*v.X + v.Y*v.Y); //length of v
		Dipole.SetLength(l);

		fePolyline<double>::V2 p = first;
		
		double phi = gCalcPhi(v.X, v.Y);
		double sinphi = 0;
		double cosphi = 0;
		FastMath::SinCos(phi, sinphi, cosphi);

		fePolyline<double>::V2 center;
		center.X = p.X + 0.5*v.X;
		center.Y = p.Y + 0.5*v.Y;
		
		for(int triind=0; triind<RecTris.GetTriNum(); triind++)
		{
			//get tri
			fePolyline<double>::V2 tri[3];
			RecTris.GetTri(triind, tri[0], tri[1], tri[2]);				
			double square = RecTris.GetArea(triind);

			//integrate eds
			fePolyline<double>::V2 p;
			p.X = (tri[0].X + tri[1].X + tri[2].X) / 3;
			p.Y = (tri[0].Y + tri[1].Y + tri[2].Y) / 3;

			double x = p.X - center.X;
			double y = p.Y - center.Y;
			double x1 = x*cosphi + y*sinphi;
			double y1 = -x*sinphi + y*cosphi;


			Dipole.DBzDt(x1, y1, iZ, DbzDt, iTimes, iTimesCount);
	
			for(int j=0; j<iTimesCount; j++)
			{
				oResV[j] += square * DbzDt[j];
			}
			
		}
	}

	return true;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
LineEdsCalcTask::~LineEdsCalcTask()
{

}

